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Abstract — The geometric heights obtained from GNSS observations cannot be used for engineering works as 
they are not reduced to the geoid. This study presents practical local geoid modelling from gravimetric 
observations using the modified Stokes integral for engineering applications in Benin City. A total of 52 points 
were observed with GNSS receivers and a gravimeter to respectively obtain their positions and absolute gravity 
values. The theoretical gravity values of the points were computed on the Clarke 1880 ellipsoid to obtain their 
local gravity anomalies. The modified Stokes integral was applied to compute the geoid heights of the points. 
The combined topographic effect was applied to the computed geoid heights of the points to obtain their precise 
geoid heights. The mean of the precise geoid heights of the points was computed to obtain the local gravimetric 
geoid model of the study area. The determined geoid model was validated for its reliability as well as the 
accuracy using the RMSE index. It is recommended that the use of assumed, as well as handheld GPS receiver 
heights for engineering works should be totally abolished as this study has established the local geoid model of 
Benin City. 
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The geoid is an equipotential (level) surface of the 
earth's gravity field which coincides with mean sea level 
(MSL) in the open oceans. As such, the geoid provides a 
meaningful reference frame for defining heights. The 
importance of accurately modelling the geoid has 
increased in recent years with the advent of satellite 
positioning systems such as the Global Navigation 
Satellite System (GNSS). GNSS provides height 
information relative to a best-fitting earth ellipsoid rather 
than the geoid (Seager et al ., 1999; Yilmaz and Arslan, 
2006). To convert ellipsoidal heights derived from GNSS 
to conventional (and meaningful) orthometric heights, the 
relationship between the geoid and the ellipsoid must be 


I. INTRODUCTION 


known (Kotsakis and Sideris, 1999; Yilmaz and Arslan, 
2006). The fundamental relationship that ties ellipsoidal 
heights obtained from Global Navigation Satellite System 
(GNSS) measurements and heights with respect to a 
vertical datum established using spirit levelling and gravity 
data, is to the first approximation given by (Heiskanen and 
Moritz, 1967; Krynski and Lyszkowicz, 2006; Oluyori et 
al., 2018) as: 


h = H +N 


( 1 ) 


Where h is the ellipsoidal height, H is the orthometric 
height, and N is the geoid undulation. Figure 1 shows the 
relationship between the orthometric, geoid and ellipsoidal 
heights. 



Fig: 1: Relationship between Orthometric, Geoid and Ellipsoidal Heights 
Source: Eteje et al. (2018) 
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Geoid Modelling is carried out by various methods such 
as the gravimetric, geometric, Astro-geodetic, 
transformation and gravimetric-geometric methods. The 
gravimetric method can be carried out by the application 
of the well-known Stokes-integral, equation (2) and the 
use of accurately determined absolute gravity data 
(Heiskanen and Moritz, 1967, Eteje, 2015 and Eteje et al., 
2018). 

V = f jj AgS^da (2) 

• <7 

Where N is geoid undulation, A g is gravity anomaly, 

S{y/) is Stokes function, y is normal gravity on the 

reference ellipsoid and R is mean radius of the earth. 

The application of the Stokes integral for local geoid 
modelling requires the solving/integration of the modified 
Stokes integral given in equation (2). Here, the geoid 
height of each point in a study area is computed and 
corrected for the combined topographic effect to obtain a 
precise geoid height of the point. The geoid model of the 
area of study is the average, as well as the mean of the 
geoid heights of the selected/observed points. The 
accuracy of the model is obtained using the Root Mean 
Square Error (RMSE) index. The RMSE of the model is 
computed by comparing the results obtained from 
physical measurements carried out on some selected 
points with their corresponding results from the 
determined geoid model. 

The geometric heights obtained from GNSS observations 
cannot be used for engineering works such as the 
installation of communication mast, road construction, 
building construction, the building of crude oil facilities, 


line and underground drainage construction, bridges, 
viaduct and tunnel constructions. The geometric, as well 
as ellipsoidal heights can be converted to practical heights 
reduced to the geoid, as well as the mean sea level. The 
conversion of the geometric heights of points to 
orthometric heights requires the knowledge of the geoid- 
ellipsoid separations of the points. Orthometric heights 
are required in engineering constructions to decide on 
suitable gradients to direct the flow of water and to decide 
on various levels of proposed constructions. With 
orthometric height, water can flow from one point to 
another while with ellipsoidal height, it is impossible. 
Benin City is the capital of Edo State. It is a developing 
area which requires a geoid model for meaningful and 
physical developments. Therefore, this study determines 
the local geoid model of Benin City from gravimetric 
observations using the modified Stokes' integral. 

1.1 The Study Area 

Benin City is the capital of Edo State in Southern 
Nigeria. It is a City approximately 40 kilometres north of 
the Benin River. The City is also linked by roads to 
Asaba, Sapele, Siluko, Okene, and Ubiaja and is served 
by air and the Niger River delta ports of Koko and Sapele. 
The City is made up of three Local Government Areas, 
Oredo LGA, Ikpoba Okha LGA and Egor LGA. It has a 
total population of 1,749,316 according to 2019 NPC 
projection. It covers a total area of about 1,204 km 2 . 
Benin City is bounded by UTM zone 31 coordinates 
660000mN and 712500mN, and 770000mE and 
815000mE. Figures 2a and b show the maps of the study 
area. 



Fig. 2 a: Map of Edo State 



Source: Ministry of Lands and Surveys, Benin City 
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1.2 Integration of Stokes’s Formula 

According to Eteje et al. (2018), using the modified 
Stokes integral given in equation (2), the geoid heights of 
points can be computed if their gravity anomalies and 


geographic coordinates are known. Featherstone and 
Olliver (1997) gave the integration of equation (2), as 
well as the Stokes integral as 


N = 


rAg 


8 y 


■6 sin 2 y/ a ln< sin 


f V a ^ 




■24 sin 3 


Vo 


+ sin" 


r \ 


vA 

2 J 


> + 16sin 


r vA 

V 2 J 


+ 12siiT 


Vo 


\\ 


\ ^ J 


12sin^ 


j 


Vo 


Acosy/ 0 + 5cos2^ o -1 


\ ^ j 


Where N is the geoidal height of individual point, y/ 0 is 
the surface spherical radius, y is the theoretical as well 
as normal gravity, A g is the gravity anomaly and r = R 


(3) 


Cp = Mean latitude of the points 
(p l = Latitude of individual point 
A = Mean longitude of the points 


is the mean radius of the earth. So, the computation of the 
geoid heights of points using equation (3) requires the use 

of the surface spherical radius, y/ 0 , theoretical, as well as 
normal gravity values, y , gravity anomalies, A g , and 
the mean radius of the earth, r = R of the points. 

1.3 Surface Spherical Radius Computation 

The surface spherical radius, y/ 0 is computed as 
(Shrivastava et al., 2015) 

cosy/ - sin^sin^ 1 + cos^cos^ 1 cosf/i 1 - A) (4) 


A 1 = Longitude of individual point 

1.5 Theoretical Gravity Computation 

To obtain the local gravity anomalies of points in a 
study area, the normal, as well as the latitude gravity, is 
computed on a specified ellipsoid. That is, the ellipsoid 
adopted for geodetic computation in the area or region of 
study. Eteje et al. (2018) gave the model for the 
computation of the theoretical gravity on the Clarke 1880 
ellipsoid as 


Where, 


STClarkel 880 ( B) 


9.78051938 


l + 0.0052474(36sin 2 </ 
-0.0000087S85sin 2 2 <p J 


ms 


(5) 


Where, 

STciarkei 880 = Theoretical gravity on the Clarke 
1880 ellipsoid 
(p = Station latitude 

1.4 Gravity Anomaly Computation 

The gravity anomaly, A g, is the difference between 
the observed gravity value (g) reduced to the geoid, and a 
normal, or theoretical, computed gravity value ( y 0 ) at the 

mean earth ellipsoid, where, the actual gravity potential on 
the geoid equal the normal gravity potential at the 
ellipsoid, at the projection of the same terrain point on the 
geoid and the ellipsoid respectively, that is (Dawod, 1998 
and Eteje et al , 2019) 

kg = g~Yo ( 6 ) 

Considering the nature of the topography of the earth 
surface, which is irregular in shape, there are two basic 
types of gravity anomalies (free air and Bouguer 
anomalies). In this study, it was only the free air correction 
that was applied. 


1.6 Free Air Correction 

This is the first step for reducing topography effects. It 
simply corrects for the change in the elevation of the 
gravity meter, considering only air (hence a free-air) being 
between the meter and selected datum. According to Aziz 
et al. (2010), this correction is added to the observed 
gravity because the increased radial distance of the station 
from the centre of the Earth results in a lower observed 
gravity value than if the station were at the local datum. 
The formula to calculate the magnitude of the reduction in 
practice is given by Eteje et al. (2019) as 

g FA =-^-H s =-308.6HjuGal 

r (7) 

= - 0 . 3086 HmGal 

Where, 

H = Station orthometric height in metres 
g = Mean value of gravity (980500 mGal) 
r = Mean radius of the Earth 

1.7 Mean Radius of the Earth Computation 

The mean radius of the earth, r = R was computed 
using: 

R = y/MN (8) 
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Where M is the radius of curvature along the meridian 
section and N is the radius of curvature in prime vertical. 
The formula for computation of the radius of curvature in 
prime vertical, N is given as (Ono, 2009) 


N = 


a 

(l-(2/-/ 2 )sin 2 ^ 


(9) 


while that for computation of the radius of curvature in 
meridian section, M is given as (Kotsakis, 2008) 


M = 


a( 1-e 2 ) 

(l — e 2 sin 2 (p)^ 


( 10 ) 


Where, 

a — Semi-major axis 

(p = Latitude of an observation point 


e 1 =2f — f 2 = Eccentricity squared (Eteje et 
al., 2019) 

r a-b 

J =-= Flattening 

a 


b — Semi-minor axis 

1.8 Combined Topographic Effect Computation 

To obtain a precise geoid height of a point, the 
combined topographic effect is calculated and applied to 
the computed geoid height of the point. The formula for 
the computation of the combined topographic effect, 

sn t c ° 0 z is given as (Sjoberg, 2000 and Kuczynska- 
Siehien et al, 2016): 

H 2 + — H 2 1 ( 11 ) 

3 R J 

Where G is the earth gravitational constant, p is density, 

R is the mean radius of the earth and H is the orthometric 
height of observation point which can be obtained from the 
DTM of the area. 

1.9 The Geoid Model 

The final gravimetric geoid model is the mean of the 
geoid heights computed with equations (3) and (11) and it 
is obtained using 

1 n 

Gravimetric Geoid Model = — ^ N l (12) 

n i=1 


xkj To p° _ 2jiGp 

Uiy Comb ~ 

r 


Where, 


N t = Geoid height of points computed using 


equations (3) and (11), 


n = Total number of points 

1.10 Accuracy of the Gravimetric Geoid Model 

The accuracy of the determined local gravimetric geoid 
model is obtained using the Root Mean Square Error, 
RMSE index. To evaluate the determined local gravimetric 
geoid model accuracy, the orthometric heights computed 
from the differences between the model geoid heights and 
ellipsoidal heights of some selected points are compared 
with their (the points) respective orthometric heights 
obtained from spirit levelling to get the residuals. The 
computed residuals and the total number of points are used 
to calculate the RMSE of the model. The Root Mean 
Square Error, RMSE index for the computation of 
gravimetric geoid model accuracy as given by Kao et al. 
(2017) and Eteje and Oduyebo (2018) is 

\v T v 

RMSE = ±J- (13) 

V n 

Where, 

^ = HObserved ~ ^Model (Residual) 

Hobserved = Observed Orthometric Height of Point 
H Model = Model Orthometric Height of Point 
n = Number of Points 

II. METHODOLOGY 

The adopted methodology was divided into different 
stages of data acquisition, data processing, and results 
presentation and analysis. Figure 3 shows the adopted 
methodology flow chart. 

2.1 Data Acquisition 

A total of 52 points were used in the study. The points 
included two primary control stations (XSU 92 and XSU 
100 were respectively located in Edo College and School 
of Nursing premises). The other 50 points were selected 
along the major roads of the City (See Figure 4). Spirit 
levelling was carried out on 3 of the 50 points for 
validation purpose. GNSS observations were carried out 
using CHC 900 dual-frequency GNSS receivers to obtain 
the coordinates and ellipsoidal heights of the points. The 
observations were carried out relative to control station 
XSU 92 using the static method (See Figures 5 and 6). 
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Fig 3: Flow Chart of the Adopted Methodology 



Fig 4: Selected GPS and Gravity Points 

The selected points were observed with a gravimeter 
(SCINTREX CG-5 Autograv Gravimeter) to obtain their 
absolute gravity values. The observations were carried by 
an expert, a Geophysicist from Mountain Top University, 
Ibafo, Ogun State. The gravity observations of the points 


Fig. 6: Rover Receiver at One of the 
Selected Points (RR01) at Ring Road 
were carried out in seven different loops relative to a 
point whose absolute gravity value was known and 
located within the Benin City Airport premises (See 
Figures 7 and 8). 


Fig. 5: Base Receiver at Control 
Station XSU92 



Fig. 7: Gravimeter Set over Reference 
Station at Benin City Airport 

Spirit levelling was carried out on the 3 validation points, 
as well as test points using control station XSU100 as a 


Fig. 8: Gravimeter at One of the Selected Points 
at Bypass along Benin-Sapele Road 
point of known orthometric height to obtain their practical 
heights. The levelling of the test points was carried out in 
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two loops. The first loop started from control station through TP2 to TP3 and closed back on TP1. Figure 9 

XSU100 to TP1 and closed back on control station shows the validation points'levelling loops. 

XSU100 while the second loop levelling started from TP1 



Fig. 9: Validation Points Levelling Loops 


2.2 Data Processing 

The GNSS observations were respectively downloaded 
and processed with HcLoader and Compass Post¬ 
processing software to obtain the positions and the 
ellipsoidal heights of the points. The coordinates and the 
ellipsoidal heights of the points were processed in Minna 
datum. The gravity observations of the points were 
processed by the expert who carried out the observation to 
obtain their absolute gravity values. All the necessary 
corrections such as drift correction, etc were applied 
during the processing. The theoretical gravity values of the 
points were computed on the local (Minna) datum ellipsoid 
(Clarke 1880 ellipsoid) using the latitude coordinates of 
the points, as well as equation (5). The gravity anomalies 
of the points were computed by finding the differences 
between the absolute gravity values of the points and their 
corresponding theoretical gravity values, as well as using 
equation (6). The computation of the free air correction 
requires the application of the orthometric heights of the 
points. And these were obtained by interpolation using the 
orthometric heights and the absolute gravity values of the 
two primary control stations (XSU 100 and XSU 96). The 
orthometric heights of the points were interpolated as there 
was no Digital Terrain Model (DTM) of the study area. 
The free air correction was applied to the computed 
gravity anomalies of the points using equation (7). The 
free air and the Bouguer gravity anomalies of the points 
were computed but the free air gravity anomalies were 
used in the study. This is because the geoid heights of the 
two primary control stations obtained from their known 
orthometric and ellipsoidal heights approximated the geoid 
heights of the stations computed using the free air gravity 
anomalies, as well as equation (3). The gravimetric geoid 
heights of the points were computed with the geographic 
coordinates, free air gravity anomalies and the theoretical 


gravity of the points using equation (3). The computation 
of the gravimetric geoid heights of the points required the 
application of the spherical radius and the mean radius of 
the earth. The spherical radius and the mean radius of the 
earth were respectively computed using equations (4) and 
(8). Also, the computation of the mean radius of the earth 
required the computation of the radius of curvature in 
prime vertical and in meridian section using equations (9) 
and (10) respectively. The computed gravimetric geoid 
heights of the points using equation (3) were co-geoid 
heights. To obtain precise gravimetric geoid heights of the 
points, the combined topographic effect has to be 
computed and applied to the co-geoid heights. The 
combined topographic effect was computed using equation 
(11). The final local gravimetric geoid model of the study 
area was obtained by finding the mean of the precise geoid 
heights of the points using equation (12). The spirit 
levelling of the validation points was reduced using the 
height of instrument method, as well as collimation 
method. The RMSE, as well as the accuracy of the model, 
was computed by finding the square root of the mean of 
the differences between the model and the observed 
orthometric heights of the validation points and the two 
control stations using equations (13). 

III. RESULTS PRESENTATION AND ANALYSIS 
3.1 Analysis of the GNSS Observation Results 

The DGPS observations were carried out to obtain the 
coordinates and ellipsoidal heights of the selected points. 
The DGPS observations were processed using Compass 
post-processing software. From the processing of the 
DGPS observations results, it was seen that the processed 
observations passed both the Network Adjustment Test 
and the X-Square (Chi-square) Test. This implied that the 
normal matrix generated was a regular one and inverted 
accordingly for the calculation of residuals. 
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3.2 Analysis of the Validation Points Levelling 

Table 1 presents the closure errors/accuracy of the two 
loops of the validation points levelling. The levelling of 
the validation points was done to obtain the orthometric 
heights of the validation points. The levelling was carried 
in two loops. The first loop started from XSU100 to TP1 
and closed back on XSU100 while the second loop started 
from TP1 through TP2 to TP3 and closed back on TP1. 


From Table 1 it is seen that the closure error of the first 
loop is -0.006m while the second loop closure error is 
0.009m which were within millimetres standard. The 
results were accepted as the closing errors of the two 
loops. The high accuracy of the levelling was as a result of 
the fairly flat topography of the study area, the observer’s 
know-how and the equipment used. 


Table 1: Known and Observed Heights of the Closing Stations 


Station 

Description 

H(known) (m) 

H(observed) (m) 

AH (m) 

Loop 1 

Starting and Closing Station (XSU100) 

76.377 

76.383 

-0.006 

Loop 2 

Starting and Closing Station (TP1) 

60.912 

60.903 

0.009 


3.3 

Analysis of the Gravimetric Geoid Model 

Table 2 and Figures 10a and 10b respectively present 
the coordinates of the selected points, their corresponding 
computed local gravimetric geoid heights and the 
determined geoid model, and their surface and contour 
plots. The gravimetric geoid heights of the points were 
computed using the integration of modified Stokes' 
integral. The gravity anomalies of the points used for the 
computation of the local gravimetric geoid heights of the 
points were computed on the Clarke 1880 ellipsoid which 
is the ellipsoid adopted for local geodetic computation in 


Nigeri 

a. From Table 2, it is seen that the determined geoid 
model (mean of geoid heights) is 2.066m. This implied 
that to convert ellipsoidal heights of points to orthometric 
heights in Benin City, 2.066m will be subtracted from the 
ellipsoidal heights of the points. Also, from Figures 10a 
and 10b, it is seen that at the centre, there are depressions 
with small cliffs closed to them while at some distances 
from the centre there are small depressions which implies 
that the geoid heights of the study area vary with no 
constant value. 


Table 2: Coordinates of the Selected Points and their Respective Local Gravimetric Heights 


STATION 

Northing 

Easting 

Free Air Geoid Height, N Corrected for 
Combined Topographic Effect 

XSU92 

257998.9800 

357763.3720 

2.086 

RR01 

257885.3227 

355124.0166 

2.420 

SR01 

254586.4919 

355927.3773 

1.588 

SR02 

253034.8393 

356093.6672 

1.978 

SR04 

249754.3940 

356486.6091 

2.520 

SR05 

245976.7564 

356615.1406 

2.802 

SR06 

244918.0916 

356628.3396 

3.266 

XSU100 

252357.6855 

356143.1412 

2.098 

AR01 

257163.2838 

354191.2450 

0.685 

AR02 

256084.6701 

352774.2792 

0.720 

AR03 

253855.5374 

351456.5724 

1.436 

AR04 

253286.3364 

351007.1375 

1.439 

UU01 

267318.6942 

352896.8470 

4.658 

UU02 

265145.3515 

353468.5482 

3.498 

UU03 

262403.8368 

354173.5295 

1.981 

UU04 

260409.1199 

354602.6925 

1.276 

UU05 

259407.1043 

355613.0973 

1.346 

UU06 

258099.7270 

356379.9681 

1.489 

UU07 

257012.1709 

355964.2081 

1.329 

UU08 

256422.9868 

355521.4167 

1.263 

AD01 

260514.8753 

359958.1194 

2.986 
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AD02 

261374.7703 

361092.6917 

4.019 

AD03 

261867.2294 

361745.9231 

4.420 

AK01 

258765.7701 

355982.5939 

1.376 

AK02 

259528.7811 

356853.3277 

1.473 

AK03 

259836.4068 

358613.7581 

2.252 

AK04 

259620.2060 

360694.5908 

3.101 

AK05 

259332.1257 

362604.6963 

3.954 

MR01 

259195.6591 

355569.5117 

1.300 

MR02 

260751.5081 

356528.1658 

1.488 

MR03 

262096.1924 

357412.2545 

1.614 

MR04 

262930.8267 

360077.3193 

4.037 

MR05 

262428.2213 

361076.8116 

4.313 

SK01 

256829.9481 

356396.3673 

1.500 

SK02 

255516.1557 

357459.1723 

2.035 

SK03 

254396.4836 

358439.3812 

2.379 

EKS 

258508.0691 

354257.9420 

0.665 

SLK0 

259220.8416 

353748.2583 

0.668 

SLK01 

259894.0672 

352909.3470 

0.781 

SLK02 

261105.4062 

351776.4441 

1.326 

SLK03 

261813.3387 

350594.2641 

1.736 

SLK04 

263367.4251 

349531.4676 

2.688 

SLK05 

264774.9356 

348869.1903 

3.357 

EK01 

257862.9575 

352479.5790 

0.078 

EK02 

257209.3523 

350068.7731 

0.983 

EK03 

255709.4653 

348058.3750 

1.729 

EK04 

254327.4139 

347366.3299 

2.001 

EK05 

252877.2407 

345740.0760 

2.516 

AIRPORT 

256224.9627 

352774.5959 

0.578 


GEOID MODEL (MEAN OF GEOID HEIGHTS) = 

2.066m 



Fig 10a: Surface Plot of Gravimetric Geoid Heights of 
the Selected Points 



346000 348000 350000 352000 354000 356000 358000 360000 362000 


Fig 10b: Contour Plot of Gravimetric Geoid Heights of 
the Selected Points 


3.4 Analysis of the Accuracy/Validation of the 
Determined Geoid Model 

Table 3 and Figure 11 respectively present the 
computed RMSE of the geoid model and the plot of the 
model and the observed orthometric heights. This was 
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done to present the consistency, as well as the reliability 
of the determined geoid model. It can be seen from Table 
3 that the computed RMSE, as well as the accuracy of the 
determined gravimetric geoid model, is 0.412m. This 
implies that ellipsoidal heights can be converted to 
orthometric height with an accuracy of 41cm using the 


determined geoid model. It can also be seen from Figure 
11 that the model and the observed orthometric heights of 
the validation points/stations are identical in shape which 
also implies the high reliability, as well as consistency of 
the determined gravimetric geoid model. 


Table 3: Observed and Model Orthometric Height, and Model RMSE/Accuracy 


STATIO 

N 

DGPS 

ELLIPSOIDA 

L HEIGHT (m) 

GEOID 

MODE 

L (m) 

MODEL 

ORTHOMETRI 

C HEIGHT (m) 

OBSERVED 

ORTHOMETRI 

C HEIGHT (m) 

DIFF. BAY 

MODEL & 

KNOWN 

ORTHOMETRIC 
HEIGHTS (m) 

DIFF. 

SQUARE 

D 

XSU92 

106.668 

2.066 

104.602 

105.441 

0.839 

0.704 

XSU100 

78.399 

2.066 
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Fig. 11: Validation Points Model and Observed Orthometric Heights 


IV. CONCLUSIONS AND RECOMMENDATIONS 

1. The local geoid model of Benin City has been 
determined to be 2.066m using the gravimetric 
method of integration of modified Stokes formula. 

2. The study has shown that ellipsoidal heights can be 
converted to orthometric heights with an accuracy of 
0.412m using the determined geoid model. 

3. It is recommended that the determined geoid model 
should be applied whenever ellipsoidal heights are to 
be converted to practical, as well as orthometric 
heights in Benin City. 

4. It is also recommended that the use of assumed, as 
well as handheld GPS receiver heights should be 


totally abolished as this study has established the local 
geoid model of Benin City. 
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